Articulation de différents langages (R, JavaScript et Python) pour la géovisualisation avec Quarto

SAGEO, 2023 - Québec, Canada

Nicolas Lambert, Timothée Giraud, Matthieu Viry, Ronan Ysebaert

07 Jun 2023

Python

  • Python : un langage polyvalent, interprété et multi-paradigme

  • De plus en plus utilisé pour la data science

  • Un écosystème robuste pour différents domaines d’application scientifiques

Écosystème pour le géospatial

Écosystème pour le géospatial - vecteur

Écosystème pour le géospatial - GeoPandas

“GeoPandas is an open source project to make working with geospatial data in python easier. GeoPandas extends the datatypes used by pandas to allow spatial operations on geometric types. Geometric operations are performed by shapely. Geopandas further depends on fiona for file access and matplotlib for plotting.”

Format des GeoDataFrame

Écosystème pour le géospatial - raster

  • Rasterio :
    • lecture / écriture de raster (wrapper de haut niveau autour de GDAL)
    • données représentées sous forme d’array NumPy
    • reprojection
    • resampling
    • virtual files
    • etc.
  • Rasterstats :
    • résumer des données raster sur la base de géométries vectorielles
    • extraction de valeurs à un point précis

Écosystème pour le géospatial - Rasterio

Exemple d’utilisation

Écosystème pour le géospatial - autre

  • Binding Python pour GRASS + Intégration dans les notebooks Jupyter

Écosystème pour le géospatial - Analyse spatiale

Ressources Python Géospatial

  • Yes

  • No

  • Maybe

Interaction R ⇄ Python dans Quarto

  • Définition de variable dans un chunk R :
a <- 42
b <- list(1, 2, 3)
c <- c(12, 13, 14)
  • Récupération depuis un chunk Python :

en utilisant la variable r, un point, et le nom de la variable R à récupérer

print(r.a)
42.0
print(r.b)
[1.0, 2.0, 3.0]
print(r.c)
[12.0, 13.0, 14.0]
  • Définition de variable dans un chunk Python :
a = 42
b = [1, 2, 3]
  • Récupération depuis un chunk R :

en utilisant la variable py, un dollar $, et le nom de la variable Python à récupérer

print(py$a)
[1] 42
print(py$b)
[1] 1 2 3

Interaction R ⇄ Python dans Quarto (suite)





  • Il est possible d’échanger des types plus complexes (data.frameDataFrame pandas, tableau numpy, etc.)

Interaction R ⇄ Python dans Quarto (suite)

  • Depuis un chunk R :
df <- data.frame(
   emp_id = c(1:5),
   emp_name = c("Rick","Dan","Michelle","Ryan","Jane"),
   salary = c(623.3,515.2,611.0,729.0,843.25),
   start_date = as.Date(c("2012-01-01", "2013-09-23", "2014-11-15", "2014-05-11", "2015-03-27")),
   stringsAsFactors = FALSE
)
  • Récupération depuis un chunk Python :
df = r.df
print(df.head())
   emp_id  emp_name  salary  start_date
0       1      Rick  623.30  2012-01-01
1       2       Dan  515.20  2013-09-23
2       3  Michelle  611.00  2014-11-15
3       4      Ryan  729.00  2014-05-11
4       5      Jane  843.25  2015-03-27
  • Depuis un chunk Python :
import pandas as pd
import numpy as np

df2 = pd.DataFrame({
  "emp_id": list(range(5)),
  "emp_name": ["Rick", "Dan", "Michelle", "Ryan", "Jane"],
  "salary": [623.3, 515.2, 611.0, 729.0, 843.25],
  "start_date": pd.to_datetime(["2012-01-01", "2013-09-23", "2014-11-15", "2014-05-11", "2015-03-27"]),
})

arr = np.array([[12, 47], [34, 90], [23, 19]])
  • Utilisation depuis un chunk R :
head(py$df2)
  emp_id emp_name salary          start_date
1      0     Rick 623.30 2012-01-01 01:00:00
2      1      Dan 515.20 2013-09-23 02:00:00
3      2 Michelle 611.00 2014-11-15 01:00:00
4      3     Ryan 729.00 2014-05-11 02:00:00
5      4     Jane 843.25 2015-03-27 01:00:00
print(py$arr)
     [,1] [,2]
[1,]   12   47
[2,]   34   90
[3,]   23   19

Python, Quarto et interactivité

  • Utilisation des widgets Jupyter possible (seulement is utilisation de l’engin de rendu jupyter, pas avec knitr - i.e. l’inverse des htmlwidgets en R qui ne fonctionne que si knitr est utilisé)

  • Exemple :

Python et géospatial dans Quarto

  • Les imports nécessaires :
import numpy as np
import rasterio as rio
import pandas as pd
import geopandas as gpd

from rasterio.warp import calculate_default_transform, reproject, Resampling
from matplotlib import pyplot as plt
from rasterio.plot import show
from rasterstats import zonal_stats
  • On ouvre les jeux de données :
reg_fp = './geom/regio_l.shp'
mun_fp = './geom/munic_s.shp'
lc_fp = '../../Téléchargements/T01_PROVINCE.tif'
lc_categories_fp = '../../Téléchargements/correspondance_raster_CL_COUV.dbf'
dst_epsg_code = 6623

regio_s = gpd.read_file(reg_fp).to_crs(epsg=dst_epsg_code)
munic_s = gpd.read_file(mun_fp).to_crs(epsg=dst_epsg_code)

Python et géospatial dans Quarto

  • Séléction de la zone qu’on souhaite étudier, en utilisant une variable de l’environnement R :
# sel = r.sel
sel = 'Québec'
  • On extrait les entités correspondantes et on les affiche, relativement au reste de la province :
extract = munic_s[munic_s.MUS_NM_MUN == sel]

ax = munic_s.plot(color="lightblue", edgecolor="grey", alpha=0.5)
ax = regio_s.plot(ax=ax, color="orange", alpha=0.8)
ax = extract.plot(ax=ax, color="red")
ax

Python et géospatial dans Quarto

  • Ouverture et reprojection d’un jeu de données raster
dst_crs = f'EPSG:{dst_epsg_code}'

with rio.open(lc_fp) as src:
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height, *src.bounds)
    kwargs = src.meta.copy()
    kwargs.update({
        'crs': dst_crs,
        'transform': transform,
        'width': width,
        'height': height
    })
    
    # Read first band
    band = src.read(1)

    # Create the destination array
    img = np.zeros_like(band)

    # Reproject and store the result in the
    # destination array
    reproject(
        band,
        img,
        src_transform=src.transform,
        src_crs=src.crs,
        dst_transform=transform,
        dst_crs=dst_crs,
        resampling=Resampling.nearest,
    )

    # Store the extent of our raster data (we will need it later for plotting)
    left = transform[2]
    top = transform[5]
    right = left + transform[0] * width
    bottom = top + transform[4] * height
    extent = [left, right, bottom, top]
(array([[255, 255, 255, ..., 255, 255, 255],
       [255, 255, 255, ..., 255, 255, 255],
       [255, 255, 255, ..., 255, 255, 255],
       ...,
       [255, 255, 255, ..., 255, 255, 255],
       [255, 255, 255, ..., 255, 255, 255],
       [255, 255, 255, ..., 255, 255, 255]], dtype=uint8), Affine(50.00000000000006, 0.0, -830584.8693000015,
       0.0, -50.00000000000006, 1303502.274700009))

Python et géospatial dans Quarto

# Replace '255' values by '0'
img[img == 255] = 0

# Create an empty figure
fig, ax = plt.subplots(figsize=(16, 12))
# Plot the raster using the previously computed extent
show(img, ax=ax, extent=extent, cmap="Set3")
# Plot Québec municipality on top
extract.plot(ax=ax, color='red', edgecolor='red', linewidth=2)

Python et géospatial dans Quarto

# Open the DBF file that contains the correspondences between the codes and the names of land use categories:
categories = pd.DataFrame(gpd.read_file(lc_categories_fp, encoding='utf-8')[['ID', 'CL_COUV', 'Descriptio']])
categories.head(11)
      ID CL_COUV                             Descriptio
0    1.0  010000                 Surfaces artificielles
1    2.0  020000                       Terres agricoles
2    3.0  060000             Milieux humides forestiers
3    4.0  070000  Milieux humides herbacés ou arbustifs
4    5.0  100000        Plans et cours d’eau intérieure
5    6.0  110101    Forêts de conifères à couvert fermé
6    7.0  110102     Forêts de feuillus à couvert fermé
7    8.0  110103          Forêts mixtes à couvert fermé
8    9.0  110200                Forêts à couvert ouvert
9   10.0  000000                         Pas de données
10  11.0  000001               En attente de traitement

Python et géospatial dans Quarto

  • Compute the zonal statistics (how many cell of each type in the selection):
stats = zonal_stats(extract, img, affine=transform, categorical=True, nodata=0)
stats
[{1: 85316, 2: 12475, 3: 2737, 4: 756, 5: 3494, 6: 3556, 7: 38023, 8: 38147, 9: 69, 10: 43}, {8: 4}]
  • Lets group the values into a single dict:
# We convert the dict objects into Counter objects, which allows us to add them up by simply using the "+" operator.
from collections import Counter
stats = Counter(stats[0]) + Counter(stats[1])
stats = dict(stats) # Convert back to dict
stats
{1: 85316, 2: 12475, 3: 2737, 4: 756, 5: 3494, 6: 3556, 7: 38023, 8: 38151, 9: 69, 10: 43}

Python

  • Lets plot an histogram of this result, using the appropriate names of land use categories:
# Use the category names from the `categories` DataFrame
x = [
  categories[categories.ID == int(v)]['Descriptio'].iloc[0]
  for v in stats.keys()
]
# and the values we just computed with the `zonal_stats` function
height = stats.values()

fig, ax = plt.subplots(figsize=(8, 5))
bar_container = ax.bar(x, height)
ax.set_xticklabels(
    ax.get_xticklabels(),
    rotation=45,
    horizontalalignment='right'
)
ax

Python